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ABSTRACT 


This report provides a general theory of the diffraction of 
ultrasonic waves. The theory is then used to find the apparent 
attenuation which would result if certain damage states (transverse 
cracks and delaminations) are introduced into a graphite/epoxy 
laminate through which the ultrasonic wave passes. 

Some experimental data for two different laminates ([0, 90, ±45]s 
and [0, ±45, 90]s) is presented which shows changes in the apparent 
attenuation of about 1 dB. These changes generally occur at loads 
which correspond to the range predicted for the formation of the 
above mentioned damage. Though no exact correlation between 
theoretical and experimental results is given, the predicted changes 
in the attenuation for several simple and common damage states are 
well within the range of experimental values. 

It is hoped that the technique described herein can be further 
developed and used to detect the formation and growth of damage in 
composite specimens in regions not readily visible by conventional 
techniques . 
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1 . INTRODUCTION 


In a real sense, the macroscopic behavior of a material is deter- 
mined by interactions which occur on a microscopic or sub-microscopic 
level. This is especially true of the failure of a material wherein 
small local failures gradually intertwine until they produce an effect 
which is catastrophic for the material, and, quite possibly, to the 
structure containing that material. 

An understanding of the failure of composite laminates, then, 
must be preceded by an understanding of these microscopic failures; 
in general, the major types of microscopic damage are transverse cracks 
and delaminations. Usually, though, only the outermost edge of a 
specimen can be observed visually, and the need to "observe" the damage 
interior to a specimen, both during and after the formation of that 
damage, calls for new methods and techniques in the field of non- 
destructive testing and engineering. 

Several techniques, both new and old, are under investigation 
at Virginia Polytechnic Institute and State University. Among these are 
vibrothermography , acoustic emissions, acoustic emission signature 
analysis, ultrasonic attenuation measurements, and replications. This 
report investigates the possibility of one more type of ultrasonic 
testing, where the formation of damage in a composite specimen serves 
as a rudimentary diffraction grating which causes a change in the 
apparent attenuation of the specimen. 

The underlying purpose here, then, is to serve as a basis for 
future investigations. This text merely purports to show that such 
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attenuation changes could be caused by the formation of damage, not 
that they are caused by such formation. That proof, and the possibl 
application of the method herein described to measuring microscopic 
failure in composites, will require considerable additional time and 
effort. 


2. LAMINATE PLATE THEORY 


The analysis of laminated composite plates under in-plane loads is 
well documented, and will be only briefly presented here. Readers are 
referred to [1] for more details. 

The simplified theory assumes that each lamina is under a state of 
two-dimensional stress, and that all displacements are small (Kirchoff- 
Love hypotheses). The in-plane displacements are given by 


u(x,y) - u°(x,y) - z Midi 
v(x.y) = v°(x,y) - z Midi 


0 ) 


where u and v are the x and y displacements at a point in the 
laminate; u° and v° are the x and y displacements of the mid- 
plane of the laminate; z is the vertical distance from the midplane to 
the point; and w is the vertical displacement of the midplane. The 
coordinate system is shown in Fig. 1. The x-y-z system is the global 
coordinate system, and is the same for all layers in the laminate. The 
1-2-3 system is known as the material coordinate system, and varies 
from layer to layer. The angle between the 1-axis and x-axis is often 
used to refer to a particular lamina; thus, a lamina whose angular 
orientation is 45° is often called a 45° layer. Material properties are 
usually expressed in the local (material) system, and then transformed 
to the global system when necessary. 

Differentiation of Eqn. (1) to find the in-plane strains yields: 
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are called the midplane curvatures. 
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stiffness 
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As in classical plate theory, the stresses and z times the 
stresses are integrated over the thickness of the plate to obtain the 
force and moment resultants. Since u°, v°, and w are independent 
of z, these integrations reduce to the integrations of the stiffness 
matrix across the thickness. The force and moment resultants, in 
terms of the strains and curvatures are 



( 8 ) 


where ^ 

< A ir B ir D ij> • ]0. z . z 2 ) Q f;j * 0) 

' ' -H 
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when H is the laminate thickness. If the laminate is symmetric about 
the midplane, it can be shown that the B matrix is identically zero, 
reducing Eqns. (7) and (8) to 
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and 



Once the midplane strains and curvatures are determined in terms 
of the applied loads, from either Eqns. (7) and (8), or Eqns. (10) and 
(11), they are substituted into Eqn. (4) to determine the stresses in 
each layer. For the case of a symmetric laminate, the stresses in 

XL 

the rn layer when no external moments are applied, are 



[Q] m [A] 


-1 




(12) 


where 
of the 




is used to denote the transformed reduced stiffness matrix 
layer. 


2.1 Curing Stresses 


Because the composite laminate is used at temperatures much lower 
than the temperature at which it was cured, residual internal stresses 
result from the different expansions and contractions of each layer. 
These must also be included in the stress analysis of the composite. 

This becomes especially important when considering damage development in 
composites, as will be discussed later. Though what follows can be 
generalized to apply to any laminate, it is developed here only for 
symmetric laminates. . 
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Since the laminate is symmetric, there will be no bending; hence, 
the strain in each layer will be the same. This strain is composed of 
two parts; the free thermal strain, and the strain imposed by the 
constraints of the other layers. Thus 

(e°> = (e T } k + {e m > k (13) 

where {e 0 } are the three midplane strains, {E T > k > the three thermal 

XL m 

strains in the k layer, and {e > k , the three mechanical strains. 
Because the mechanical strains arise from the constraints of the other 
layers, they determine the stresses in each layer, thus 

[{e 0 } - te T > k ] 04) 

The thermal strains have the form of 



where aT is the temperature difference causing the thermal strains, 

• 1L 

and the a's are the coefficients of thermal expansion of the k ir 
layer. 

If there is no net force acting on the laminate, the integral of 
the stresses across the laminate thickness must be zero . Again, noting 
that {e 0 } is independent of the z direction, this integral becomes 
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H 

2 

' [ID k {a> k AT dz = 0 


(16) 


-H -H 

2 2 

Calling the last integral of Eqn. (16) {N , Eqn. (16) is solved 
for {e 0 } to give 

{e°} = [A]' 1 W T (17) 

Equation (14) is then used to solve for the curing stresses in each 
layer, 

2.2 The Pipes-Pagano Interlaminar Normal 
Stress Approximation 

Laminate theory assumes that the entire plate is under a state 
of two-dimensional stress, and this is essentially true in regions 
away from the free edges a distance approximately equal to the thickness 
of the laminate. In these boundary regions, however, significant 

interlaminar normal stresses ( a ) exist. Pipes and Pagano [2] have 

z 

proposed an approximate solution for the normal stress in the boundary 
region. No details are given here, but the interested reader is invited 
to refer to their paper. The greatest importance of their approximation 
is the relative magnitude of the normal stress and its sign (tensile or 
compressive) through the thickness of the laminate. This will be 
discussed in more detail later. 


3. DAMAGE DEVELOPMENT IN COMPOSITE MATERIALS 


The subject of damage development in composites has received much 
attention, but the phenomenological process remains unexplained for the 
most part. The discussion presented here is largely conjecture, although 
supported by some experimental observations. 

Two different composite laminates are discussed here. The first, 
hereafter called type I, has a [0°, ±45°, 90°] s stacking sequence. The 
second, hereafter called type II, has a [0°, 90°, ±45°] s stacking 
sequence. Average lamina properties of the AS-3501 graphite epoxy used 
in constructing the specimens are presented in Table I. 

The two most frequently observed types of damage in composites are 
transverse cracks and delaminations. Transverse cracks are those which 
exist parallel to the fiber (1) direction, caused mainly by stresses 
in the 2 direction (a^). Delaminations are separations between laminae 
caused by interlaminar normal stresses (ag). Delaminations are generally 
confined to a region close to the free edges of a specimen, while 
transverse cracks are thought to extend across the width of the specimen. 
Though seemingly unrelated, both types of damage occur because of the 
relatively low strength of the epoxy used to hold the graphite fibers 
together. 

A computer program was written to calculate the curing stresses 
for a symmetric laminate, as well as to predict the average laminate 
properties and the applied resultant (N x ) necessary for the failure of at 
least one lamina. The failure criterion used was a Tsai -Hill criterion 
[1] which states that failure occurs when the following condition is 
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met: 
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0 1°2 . °2 . 
? ~ 


-J2 = ! 

o o * 


(18) 


The values used for X, Y, and S are also presented in Table I. The 
curing, mechanical, and total stresses at first-ply failure are 
presented in Table II. Note that the stresses and average laminate 
properties listed are for both types of laminates, since laminate theory 
cannot distinguish between their in-plane behavior. Table III presents 
the results of the Pipes-Pagano approximation of the maximum interlaminar 
normal stress existing between laminae. Note that these stresses are 
different for the two laminates. The applied load resultant, N , is 

A 

875 Ibs/in for both tables. 

For both laminates, the 90° plies fail first, and failure occurs 

y 

when ct 2 is approximately equal to Y. The curing stress (o^) for 
these laminae is 4.77 ksi, and if N is 875 Ib/in, the total stress 

X 

(a 2 ) in the 90° layers is 8.2 ksi, resulting in a predicted total 
failure of the 90° plies by transverse cracking. Experimentally, 
cracks have been observed to begin at about 900-1000 Ibs/in; however, 
the entire ply does not fail at once, but continues to develop cracks 
at loads up to about 2000 Ibs/in. This gradual failure can be explained 
in part >by a spatial variation of the value of Y, since the 8.2 ksi 
used is only an average value. A better explanation, though, is 
contained in a paper by Pagano [4], where he reports that the absorption 
of around 0,3 % moisture, by weight, by the composite will almost 
completely relax the thermal stresses. If these stresses are relaxed. 
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TABLE I 

SUMMARY OF PROPERTIES* FOR AS- 3501 


Elastic 


Sti ffness 
Properties 

Tensile 

Value 

Compressive 

Value 

E i 

20.2 x 10 6 psi 

17 x 10® psi 

E 2 

1 .4 x 10® psi 

1.6 x TO® psi 

v 12 

0.28 

0.28 

G 12 

0.65 x 10 6 psi 

0.65 x 10® psi 

Strength properties 

X 

235,000 psi 

180,000 psi 

Y 

8,200 psi 

25,000 psi 

S 

17,900 psi 

17,900 psi 


Other properties 


Thermal expansion coefficients 


-6 


a-| (fiber direction) = -0.2 x 10 in/in/ F 


a 2 (transverse direction) = 13 x 10 
Fiber volume fraction 0.62 
Ply thickness 0.0052 ± 0.0004 in 
Void content 2 % 

Density 0.057 lb/in'* 

Stress-free temperature 278°F 

*Data supplied by Hercules, Inc. 


-6 , 


in/in/ F 
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TABLE II 

THERMAL, MECHANICAL, AND TOTAL STRESSES, AND AVERAGE 
LAMINATE PROPERTIES AT FIRST PLY FAILURE IN TYPE I AND II LAMINATES 


Layer 

Orientation 

Stress 

Direction 

Thermal 
Stress(ksi ) 

Mech. 

Stress(ksi ) 

Total 
Stress (ks 

0° 

a l 

-4.77 

53.9 

49.12 


a 2 

4.77 

- 0.07 

4.70 


t 12 

0 

0 

0 

90° 

a l 

-4.77 

-15.17 

-19.94 


a 2 

4.77 

3.42 

8.19 


T 12 

0 

0 

0 

45° 

G 1 

-4.77 

19.36 

14.59 


a 2 

4.77 

1.67 

6.44 


t 12 

0 

- 2.90 

- 2.90 

o 

LC> 

1 

a l 

-4.77 

19.36 

14.59 


a 2 

4.77 

1.67 

6.44 


x 12 

0 

2.90 

2.90 


E 

= 8.11 x 10 6 

psi 



E = 8.11 x I0 6 psi 

yy K 

V = 0.299 

xy 

G yv/ = 3.12 x 10 6 psi 
.. a y 
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TABLE I IT 


THERMAL, MECHANICAL, AND TOTAL INTERLAMINAR NORMAL 
STRESS AT FIRST PLY FAILURE IN TYPE I AND II LAMINATES 


Thermal Mech. 

Interface Stress(ksi ) Stress (ksi ) 


Total 

Stress(ksi ) 
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the applied load necessary to cause failure of the 90° layers increases 
to approximately 2000 lbs/in. 

Another observation not explained by laminate theory, or the 
application of a failure criterion to the stresses predicted by 
laminate theory, is the almost regular spacing of the transverse cracks. 
Reifsnider [3] has addressed this problem by assuming that, in the 
region of the transverse cracks, the surrounding layers take up the 
load previously carried by the cracked lamina, and attempts to find at 
what distance away from the crack the stress in the failed layer 
recovers its original value. He then assumes that another crack will 
form at this distance, and thus predicts the crack spacing. For a 
type I laminate, he predicts a spacing of 0.029" in the 90° layers, 

II H 

with observed spacings of 0.024 - 0.059 . For a type II laminate, 

o " 

he predicts a spacing in the 90 plies of 0.014 , with observed spacings 

II II 

ranging from 0.009 to 0.017 . The observed spacings depend upon both 
the maximum load and the total number of cycles that the specimen 
undergoes, and decreases with both increasing load and increasing 
number of cycles. 

One last observation worth noting is that the 90° laminae in the 
type I material act, essentially, as one layer of twice the thickness. 
Hence, the transverse cracks in type I laminates would have a much 
wider opening than those in type II laminates, where the 90° laminae 
are independent. 

Once the 90° laminae begin to fail , laminate theory is no longer 
applicable, and attempts to predict the failure of the 45° laminae 
are somewhat dubious. Experimentally, the 45° plies adjoining the 90° 
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plies begin to develop transverse cracks shortly after their development 
in the 90° layers, and the damage progresses throughout the interior 
of the laminate. 

The existence of the second type of damage is predicted by the 
Pipes-Pagano approximation. Due to the stacking sequence, very 
high tensile stresses are present between the layers of the type I 
laminate, while the stresses are predominantly compressive in the type 
II laminates. This is evidenced by extensive del ami nations along the 
edges of type I specimens. Occasionally, the type II specimens 
delaminate near the end tabs, but this appears to be a grip related 
phenomenon. Little work has been done to date to determine the depth 
of the delamination as a function of load or number of fatigue cycles, 
but it appears to extend at least a quarter of an inch into the width, 
at or near specimen fracture, of a 1" wide specimen. 

As specimens of either material are loaded in simple tension, the 
stress-strain curve becomes bilinear at loads of approximately 2000 
Ibs/in. This knee, though not evident in all specimens, is probably 
related to the amount of transverse cracking in the various layers. The 

V 0 

predicted change in E with complete failure of the 90 layers is 

XX 

about 4 %, while if all the interior layers fail, the change is about 
22 %. The actual values found range from about 6 % to 18 %, indicating 
that the knee is some combination of failure in the 90° plies as well 
as in the 45° plies. The knee in type II laminates is generally more 
gradual than in type I, for unexplained reasons. 
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4. WAVES IN SOLIDS 


Wave propagation in all materials is governed by the equation of 
motion : 



(19) 


If the material in question is linear elastic, homogeneous, and 
isotropic, the stresses can be replaced with spatial derivatives of the 
displacements, u i , and the Lame constants, A and p v to obtain 


(x + w) u j,jT + uu i,jj • (20) 

Generally, a wave is considered to have a sinusoidal time variation, 

and complex notation is used. It must be remembered, in what follows, 

that the actual wave is mathematically described by the real part (real 

part convention) of the complex wave function. Thus u. is represented 

- i wt 

by a time varying portion, e , multiplied by a spatially varying 
portion, u\ . The time portion of Eqn. (20) is then eliminated, and 
the remaining partial differential equation is solved by appropriate 
means. 

Three of the most commonly discussed waves which can propagate 
in a solid are the longitudinal plane and spherical waves, and a 
transverse plane wave. A mathematical description of any plane wave is 

X u, - A ie i(Kn A - wt > : (21) 

where n are the components of the unit vector normal to the wave 
front; w is the angular frequency of the wave; and k is the wave 
number. Numerically, k is equal to the rati o of the angular 
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frequency to the appropriate wave speed, while the longitudinal and 
transverse wave speeds, and respectively, are given by 

1/2 


v L * (^4 
. /U\V2 


( 22 ) 


V T ■ (p 


If Eqn. (20) is solved in spherical coordinates, with u^ being 
the only non-zero displacement and all quantities being independent of 
angular orientation, then a longitudinal spherical wave can be found 
which propagates with displacements given by 



Ak ( 


kr 


.2 2 

k r 


) e 


i ( kr - wt) 


(23) 


where r is the distance to a field point from the source of the 
spherical wave, and k is the longitudinal wave number. For both 
kinds of longitudinal waves, the particle displacements are always 
parallel to the wave front normal, while for transverse waves, the 
displacements are perpendicular to the wave front normal. 

Equations (21) and (23) are mathematical descriptions of waves in 
ideal materials, but all real materials absorb some of the energy 
carried by the wave as it travels through the medium. This absorption 
process is called attenuation, and is usually described mathematically 
by allowing the wave number to become complex; i .e. , 

k = k + ta,"^ (24 v 


with k being the usual wave number and a the amplitude attenuation 
coefficient. Eqns. (21) and (23) then become 
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A.e 

J 


■ctfrx,, 
£ l 


gi ( kn £ x p - wt ) 


(25) 


u r = Ake-“ r e’< kr - (^ + 4?) (26) 

k r 

The attenuation coefficient is a function of many things, among which 
are: the material serving as a sound medium, the frequency w, and the 

damage or material discontinuities existant in the material. Hence, 
it is generally possible to obtain a measure of the amount of damage 
introduced into a material by external forces by measuring the change 
in attenuation. 

This coefficient is usually measured by a technique called the 
pulse-echo method. An ultrasonic transducer is fastened to one surface 
of a block of material by some bonding agent, usually a thin grease 
such as glycerin. The transducer is electrically coupled to a pulse 
generator-receiver which emits high energy bursts, or pulses, of very 
high frequency (2 - 50 MHz) and short duration (2 - 5 ysec). These 
electrical pulses are converted to mechanical vibrations (sound) of the 
same frequency by the transducer. The sound pulses then travel through 
the material and reflect off of the opposite surface, returning to the 
transducer and causing it to vibrate. A portion of the mechanical 
energy of the transducer is converted back into electrical energy 
which passes to the receiver and is displayed upon a CRT as an echo. 

The reflections continue until they eventually die out, resulting in 
a series of echoes on the CRT. An examination of the maximum 
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amplitude in each echo will yield the attenuation coefficient. 

Suppose that the block being tested is of length L. Then the 
distance traveled by each echo is 2nL, where n is the number of the 
echo. Suppose also that the wave travelling in the material can be 
adequately described by Eqn. (25). The maximum amplitude of the n tn 
echo is then given by 

u n = V 2 " aL (27) 

with A q being the initial amplitude of the sound wave. If the ratio 
of u m to u n is taken, with n being greater than m, one obtains 


u m _ e 2a(n-m)L 
u n ’ 


Hence, the attenuation coefficient is found to be 


. a = 


1 

2L(n-m) u 


u 

In - m 

n 

n 


(28) 


(29) 


with units of nepers/cm. 

The length of the specimen is important because the pulse must 
be shorter than the time required for an echo to return. If not, the 
returning echoes interfere with the original pulse, as well as with 
each other, making the determination of the maximum amplitude either 
incorrect or impossible. As discussed by Truell [5], other factors 
can also affect the accuracy of attenuation measurements. Foremost 
among these are losses which occur because of the coupling agent, 
losses resulting from non-parallelness of the two reflecting surfaces 
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of the specimen, and losses associated with the diffraction of the 
sound wave as it leaves the circular transducer. Some, or most, of 
these errors can be accounted for, as will be discussed later. 

Units other than nepers/cm. are often used to describe attenuation. 
Two of the most common are decibels per inch (dB/in) and decibels per 
echo (dB/echo). They are defined as 


20 m 

“ 2L(n-m) lo9 10 u„ 


(30) 


and 
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log 


m 


10 u 


(31) 


Unless otherwise specified, units of dB/echo will be used in what 
follows. 

It is unfortunate that most composite specimens are not very 
thick; the usual eight-ply laminate has a time between echoes of about 
1 ysec, while the minimum pulse width achievable with presently 
existing equipment at VPI & SU is between 1-1/2 to 2 ysec. For this 
reason, a fused silica delay block is used to obtain a larger time 
separation between echoes. 

Several reports ([6], [7]) have been made concerning the use of 
buffer blocks in making attenuation measurements. Papadakis 1 technique 
[6], though very useful on thin specimens, places severe restrictions 
on equipment performance which the presently used system is not able 
to meet; one in particular is the requirement of a very fast-rising 
pulse. Reference [7] relaxes these requirements, but in order to do so, 
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places a lower limit on the number of plies in the laminate at about 
fifteen. Both techniques placed the buffer block between the specimen 
and transducer, and kept track of the multiple reflections which occur 
at the buffer block-specimen interface. After thoroughly exploring the 
above techniques, they were discarded as being useless for composite 
specimens of eight plies or less, at the present time. 

Another technique which also has been tried at Virginia Tech 
consists of reversing the locations of the specimen and buffer block. 

A picture of the transducer and buffer block mounted on a specimen 
is shown in Fig. 2. The c-clamp is used to hold the three objects 
together as the specimen is being loaded in tension or fatigue. 

Figure 3 is a photograph of a typical oscilloscope trace using this 
arrangement. Note that the echo train consists of several large 
echoes (major echoes) evenly spaced, with a series of smaller echoes 
(minor echoes) trailing each major echo. These minor echoes return 
so closely to the major echoes as to be individually indistinguishable; 
hence, only the major echoes can be used to obtain any useful information. 
Figure 4 is a schematic diagram which shows the paths taken by the 
first two major echoes. The paths taken by the minor echoes have been 
omitted for clarity. These major echoes correspond to the portion 
of the wave which passes completely through the specimen and into the 
buffer block, after which it simply reflects back and forth in the 
buffer block, releasing a portion of its energy into the specimen at 
each reflection at the buffer block-specimen interface. It is this last 
portion of the wave which reaches the transducer and is shown on the 
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Figure 3. Typical Oscilloscope Trace for the 
Transducer - Specimen - Delay Block 
Combination 
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CRT as a major echo. 

To obtain a good estimate of the relative amplitudes of the first 
two major echoes, assume that a plane wave is emitted from the trans- 
ducer. Let a.j and a 3 be the amplitude attenuation coefficients of 
the specimen and buffer block, respectively, in units of decibels 
per round trip in either. As shown by Hayford [7], the maximum 
amplitude of a wave reflected from a bond between two materials is a 
reflection coefficient, R, times the maximum amplitude of the wave 
incident upon the interface. The transmitted portion is treated 
similarly, with its magnitude being equal to a transmission coefficient, 
T, times the incident amplitude. Since, in general, the side of the 
interface from which the wave is incident is important, both R and 
T will be given two subscripts. The first denotes the incident 
material. ' The second denotes the material on the other side of the 
interface; thus, R^ 3 denotes the reflection coefficient of the wave 
incident upon the 1-3 interface, from material 1. A schematic 
representation of the necessary reflection and transmission coefficients 
is shown in Fig. 5. From [7] 

T = 3 __ {32) 

31 (z 1 + z 2 )(z 2 + z 3 )e" 1f3 ^ + (z 2 - z 3 )(z 1 - z 2 )e n 


(z 3 - z 2 )(z 1 + z 2 )e" 1BA + (z 2 + z 3 )(z 2 - z-j)e lpJl 
(z 3 + z £ ) ( z -j + z 2 )e" 165, + (z 2 - z 3 )U 1 - z^e 13 * 


with 3 being the wave number in the bonding agent, and a the bond 
thickness. The acoustic impedance, z, of a material is given by the 
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product of the wave speed and density. The transmission and reflection 
coefficients and ^3 can be obtained by interchanging z-j and 

z 3 in Eqns. (32) and (33), respectively. Note that the coefficients 
are complex, indicating the possibility of a phase shift other than 
zero or pi, which is the case for an interface between only two 
materials. 

An analysis of the schematic representation of Fig. 4 is now 

given. In what follows, the word "amplitude" shall refer to the 

maximum amplitude of an echo. Assume that the amplitude of the wave 

just after leaving the transducer is given by A q . This wave travels 

through the specimen until it reaches the interface, at which point 

-ou/40 

it has an amplitude of A Q 10 . The portion that is transmitted 

is the only part of interest, and just on the other side of the 

-a-,/40 

interface, its amplitude is A q T -j 3 ^ 0 . The wave travels through 

the delay block until it strikes the lower end, where it completely 

reflects. When it again reaches the interface, the amplitude is 
-a-j/40-a 3 /20 

A o T-j 3 10 . Now this wave breaks into two parts. One portion 

is transmitted across the bond and travels through the specimen to the 

-a-,/40 

transducer, introducing a T^ and another 10 , Thus A^ , the 

amplitude of the first major echo, is given by 

-a r /20-a,/20 : ' 

■ ■"ywf 0 : : : : f 34 ) 

Continuing in the above manner and following the above-mentioned 
portion which reflects at the interface, it is found that A^, the 
amplitude of the second major echo, is given by 
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Fig. TO Attenuation vs. Load for Specimen 2626-10 (Type II) 
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A 2 = 


T 1 3 T 31 R 31 


10 


•a-j/ZO-a^/lO 


(35) 


If Eqn. (31) is now used to find the loss in dB/echo, one obtains 

a = - 20 log | Rgi | (36) 

Thus, as is seen in the above analysis, the only information which 
appears to be contained in the major echoes is the attenuation of the 
quartz buffer and the magnitude of the reflection coefficient, R^-j . 

Tests of this nature have been conducted on specimens undergoing 
quasi-static tension, and significant changes, which are either 
positive or negative, have been observed at load resultants ranging 
from 1000 to 2000 Ibs/in, as shown in Figs. 6 through 10, for both 
type I and type II specimens. It is interesting to note that these 
changes are not monotonic and can be either gradual or very sudden. 

Note also that the above load ranges span the values at which the 
90° layers are predicted to fail, and at which delaminations are 
observed to form. Obviously, these rises have nothing to do with 
changes in a^, since no loads are applied to the quartz delay block. 

It appears, then, from Eqn. (36), that the change in attenuation must 
occur because of changes in the bond (more specifically, the bond 
thickness), since there is no other quantity which could change during 
loading in this equation. That possibility will now be investigated. 

From Eqn. (33), the magnitude of R 3 -j is 
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-2z ] z^z 3 + cos 2 3*(z 2 + z 2 )z 2 + sin 2 0Ji(z2 + z\z\) 
2 z-jZ 2Z 3 + cos 2 3i(z 2 + z 2 )z 2 + sin^BUz^ + z 2 z 2 ) 


1/2 

(37) 


while the values of z^ , z^» and z 3 are 


z-j = 4000 g/cm 2 -sec 

z 2 * 2400 g/cnf-sec (38) 

z 3 = 13000 g/cm 2 -sec 

Figure 11 shows a plot of both |R 31 | and the bond loss term, 

-20 log I R31 [ * as a function of 35,. Note that I R31 ! has a maximum 
at 3& equal to ir/2, resulting in a minimum in the bond loss. This 
minimum is calculated to be 1.9 dB. Returning to Eqn. (36), a minimum 
value of a is actually a function of two variables, a 3 and R^ . 
However, fused silica is used as a buffer material precisely because 
of its very low intrinsic attenuation. Hence, it is reasonable to 
assume that the minimum value of a occurs at the minimum value of 
-20 log | R 33 1 . At this point, equals ir/2. Recalling the 
definition of the wave number, it is found that 


2irf£ _ it 

v b "2 



where f is the frequency of minimum attenuation, and is the 
velocity of sound in the bond. For glycerin, this latter value is 
approximately 2500 m/sec. Hence, the bond thickness, in inches, is 
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A specimen was loaded in tension to 2100 lbs., and values of the 
frequency at minimum attenuation were recorded at approximately 200 lb. 
increments. At zero load, this frequency was 25.5 MHz, and at 2000 lbs., 
it was 24.4 MHz, corresponding to an initial bond thickness of .00096 
in. and increasing to .00101 in., a change of .00005 in. If the 
frequency is reduced to the normal testing frequency of 5 MHz, the bond 
losses would then become 4.7 dB and 4.6 dB, respectively, indicating 
a decrease in attenuation of 0.1 dB over the range of the test. Of 
course, these values are only approximate for other tests, but they 
are reasonable estimates based upon the above assumptions. The important 
thing, however, is not the actual values, but the tendency of the bond 
to become thicker as the specimen thickness decreases when loaded. The 
bond thickness is always close to .001 in., and at this thickness, 
the attenuation due to bond losses will always decrease with increasing / 
load. Thus, it is hardly likely that bond losses can account for the 
large (1 to 2 dB) increases found in the tests run. 

While making attenuation measurements with high accuracy, 
experimentalists often make what could be called second-order corrections. 
These second-order Tosses, usually less than one decibel, are generally 
of a lesser magnitude than primary losses, such as the bond loss 
described above. Among these second-order losses are those which occur 
in the transducer, and those due to the diffraction exhibited by any 
finite-sized wave source, such as a transducer. The former, since 
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it is not affected by loads which occur in the specimen, will be 
ignored, and a correction term will be added to Eqn. (36) to express 
the diffraction losses. Thus 

a = a 3 - 20 log | R 3 -j | + “ d (41 ) 

where a D is the, as yet unspecified, diffraction loss. 

In the study of optical diffraction, three main types of 
diffraction effects are discussed. They are diffraction by single- 
and multiple-slit gratings, and diffraction by sharp edges. Sound 
waves will behave similarly, and it is readily apparent that both 
transverse cracks and delaminations in a specimen could serve to act 
as these types of diffraction gratings for ultrasonic waves of appropri- 
ate wavelengths. Thus, a short presentation of the theory of diffraction 
is in order, and is presented in the next section. 


5. ULTRASONIC DIFFRACTION THEORY 

In order to find a mathematical description of a sound wave after 
it has passed through some aperture, one need only apply the equations 
of elasticity, together with appropriate boundary conditions, and 
solve. Unfortunately, the mathematics involved has not seen fit to 
cooperate with the theoretician, necessitating the use of simpler 
equations. The theory of light diffraction has been well explored 
(see, for example, [8], [9]), and the equations derived for the scalar 
diffraction of light can be easily adapted to the diffraction of 
ultrasonic waves in certain media; notably, compressible, inviscid 
fluids. The term scalar diffraction is used for those cases where the 
vectorial aspect of light waves is insignificant. 

The use of a fluid medium to study diffraction effects in solids 
has been justified by several researchers. Truell [5] and Papadakis 
[10] have both done considerable research on diffraction effects in 
velocity and attenuation measurements in solids using a fluid medium 
for their theoretical calculations. Both showed excellent agreement 
between theory and experiment. Papadakis also modified the wave 
number, k, to account for the anisotropy of crystalline solids 
while still assuming a fluid medium, and achieved good agreement 
between theory and experiment. 

Wave motion in fluids is similar, in some respects, to that in 
solids, and is also governed by the equation of motion, Eqn. (19). 

But Eqn. (20) is simplified because, for inviscid fluids, the Lame 
constant u is zero. Thus Eqn . (20) becomes 
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Xu 


J.Jl 


= pu. 
1 


(42) 


By taking one more derivative with respect to the spatial coordinate 
, one obtains 


Au j,jii pu i,i 


The dilatation A is defined as u^ . . Hence 


XV A = pA 


(43) 


(44) 


If the dilatation is assumed to have a sinusoidal time variation, one 
obtains Helmholtz's time-independent wave equation 

(v 2 + k 2 )A = 0 (45) 

with k being the usual wave number, and the velocity of sound being 

V L - (f) 1/2 (46) 

Lambda is usually replaced by B, the bulk modulus, in fluids. For 
solids, B is given by 

B = X + | y (47) 

while for fluids, B and x are the same, since y is zero. 

The normal Hooke's Law for solids is modified slightly for fluids. 
Instead of relating strain to stress, the dilatation is related to 
the hydrostatic pressure, p, by the bulk modulus, to give 


p = -BA 


(48) 
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and if B is independent of position, Eqn. (45) becomes 


(v I 2 + k 2 )p = 0 


(49) 


To make the transition from solids to fluids complete, it is 
necessary to define attenuation in terms of a decrease in the dilatation, 
rather than in the displacement. The attenuation coefficient a is 
then 


a = 


1 


X 2~ X 1 


In 


1 


(50) 


in units of nepers/cm., where A-j and a 2 are the maximum dilatations 
at distances x-j and x 2 , respectively. Since the dilatation is 
proportional to the pressure, Eqn. (50) is also written 


a = 


i i i 

x 2~ x l p 2 


(51) 


The other definitions of attenuation, Eqns. (30) and (31) are adjusted 
accordingly. 

Kirchoff's diffraction integral ([7], [8]) was derived for scalar 
light waves, but can easily be adapted to sound waves in fluids. There 
are several methods for approaching this problem; the following is 
based upon the presentation of Sommerfeld [8]. 

Assume that some monochromatic sound disturbance, henceforth 
denoted by v, exists in a region of space V, By Green's theorem 

I (uv 2 v - vV 2 u)dV = | ( u " v |^)dS (52) 

V S 

as long as u and v are suitably continuous functions. S is the 
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surface surrounding V, and ~ is given by the dot product of n 
and the gradient of u, with n^ being the outward normal to S. 

All derivatives up to the second must exist for both u and v 
everywhere inside and on S. 

Since v is a wave variable, such as pressure or dilatation, 
it must satisfy Helmholtz's equation 

2 2 

V v = -k v (53) 


Suppose that u is also from the class of functions which satisfy 
Helmholtz's equation. Then the left-hand side of Eqn. (52) is 
identically zero, and one obtains 


< u f - v £>« 


3U< 


= 0 


Let u be the function 

ikr 

e 

u = 


(54) 


(55) 


where r is the distance from some point Q inside of S to any 
other point in or on S. Since u has a singularity at the point Q, 
a volume surrounding Q must be excluded from V. This is done by 
enclosing Q with a small sphere of radius e, whose surface is 
denoted by S'. The entire volume, along with the various surfaces 
and normals, is shown in Fig. 12. For this choice of u, Eqn. (54) 
becomes 


.ikr 


(H_ iV _ v 

V v» ■■ 


.ikr 


J 

S+S 


3n 


9n 


-)dS = 0 


(56) 
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or 


( e 1kr av 3 e lkl "> 

v 3^— )dS 


'I 


, a e ikr ikr 

(v lr V" V~ w )dS 


(57) 


3 3 

On S' , ~ = - — , and since r * e on S', Eqn. (57) becomes 


ar 


r Jkr k- „ Jkr f 

( < V f - v k V-> dS = 1 1 

S S' 


■e^av 

e ar 


.ike 


- v 


(1k-l)]dS (58) 


Because S' is a sphere, dS can be written as 


dS - e^dft 


(59) 


where n is the solid angle describing the sphere. Equation (58) then 
becomes 


4tt 

(^1 ' v “ I - ikv > + ve’ kc ]dsi 

0 


(60) 


Since the left-hand side of Eqn. (60) is independent of e, the limit 
can be taken as e approaches zero. Noting that v approaches Vg 
(the value of v at point Q) as e approaches zero, Eqn. (60) 
reduces to 


4ttv 




e ikr av a p ikr . 
e 3V - V hr ^r~) dS 


Q J v r an ' an r 
S 


(61) 


where S is any surface which completely encloses the point Q. 

Suppose, now, that v represents some monochromatic wave striking 
a rigid plane barrier with an aperture. For instance, the radiation 
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of a transducer can be thought of as a plane wave, radiating through 
some aperture whose shape is the same as the shape of the transducer. 

Let the plane of the aperture be one portion of S, and let the portion 
of a large sphere below the aperture, centered at Q, be the remainder 
of S. As shown by Sommerfeld [9], the integral of Eqn. (61) over 
this large sphere tends toward zero as the radius of the sphere increases 
without bound. However, that same result can be argued without 
resorting to the complex mathematics involved, as shown by Born [8]. 

His approach was as follows. 

A strictly monochromatic wave would have no beginning or ending 

in time, but a departure can be made from strict monochromaticity as 
long as the value of Vg is examined a long time after the beginning 

of the wave. Because the wave must travel from the aperture with a 
speed no greater than the velocity of sound in the fluid, there exists 
some length R, measured from point Q, beyond which the wave has 
not traveled. Let the radius of the large sphere be , where R^ 
is greater than R. Then all along that part of the sphere which makes 
up the lower portion of S, both v and are zero. Hence, the 
integral of Eqn. (61) is zero all along this boundary. Let the 
remainder of S be broken into two parts, with a denoting the 


aperture, and o' the rigid barrier. Then Eqn. (61) becomes 

ikr 


4 TTV g = 


(elLl ii - V i- 6- 

r 3n 9n r 


-)dS 


(62) 


a+0 ' 


At the barrier, it is reasonable to assume that both v and •— 
are zero;, thus 
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on o' 


(63) 


and at the aperture 

v = v' 

on a (64) 

9V _ 3V_‘ 

an ” an 

where v 1 denotes the value that the incident wave would have if the 
barrier were not present. Because of the boundary conditions expressed 
by Eqn. (63) the integral over o' is zero, and the value of v at Q 
becomes 



= 1_ f • ( e iKr av' 
4ir ■ r an 


a 




(65) 


Actually, as shown by Sommerfeld [9], the boundary conditions 
contained in Eqns. (63) and (64) are mathematically inconsistent, since 
Eqn. (63) implies that v is zero everywhere. These boundary conditions 
(Eqns. (63) and (64)) are, however, physically reasonable, even though 
they are only approximately correct. They are more accurate for high 
frequencies than low frequencies, and can be used to obtain reasonable 
estimates of v^ when the frequency is sufficiently high, 

Sommerfeld has modified the auxiliary function used in deriving 
Kirchoff's diffraction integral in order to eliminate the mathematical 
inconsistency, but his modification is itself no further justification 
of the assumed boundary conditions. That is, the inconsistency is 
removed, but the boundary conditions are still only approximately true. 



because the pressure on the downstream side of the baffle is very close 
to, but not exactly, zero. 

To do this, Sommerfeld replaces the auxiliary function u with 
the Green's function, G, which has the following properties: 

(v 2 + k 2 )G = 0 in V 

G = 0 on a and a' (66) 

G •+ u as r + 0 

Using the first and last conditions, Eqn. (62) becomes 



(G 


av 

an 


<r+a' 


v'|£)ds 

an 


(67) 


as before, and the second condition reduces the above equation to 



aG 

Tn 


dS 


0 + 0 * 


( 68 ) 


Thus, the boundary condition on t— is not needed, and the mathematical 

Oil 

inconsistency is avoided. If v is again assumed to vanish on o', 
then Vq is 



1_ 

4tt 




(69) 


The function G is easy to formulate if the aperture lies in 
a plane. Shown in Fig. 13 are the aperture, the point Q, a point P 
in the aperture, and a point R which is the image of Q. The point 
0 is the origin of the f>, n, and x, coordinate system which 
describes the aperture (c = 0). The coordinates of Q are (x, y, z). 
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Figure 13 

Relationships Defining Green's Function, G 
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and of R, (x, y, -z). The function G is then defined as 


G = 


ikr 1 ikr ? 
e e 6 


1 


(70) 


where 


r 2 = U - x) 2 + (n - 


y) 2 + (? - z) 2 


(71) 


and 


r 0 = (c - x; + (n - 


y) 2 + (; + z) 2 


(72) 


9 3 

For the surface o, — = - — , thus 


3G . 

3n 

Substituting, one obtains 


_ 90 3 r l _ 3G 9 r 2 
3r-| 9? " 3r 2 H 


ikr 


= 2ze (ik _ m 

9n 2 UK r ; 


(73) 


(74) 


with r-| = r 2 = r when c = 0. 
Hence, Vg is given by 


V Q 2tt 


.ikr 


Y~ Ok - p* s 


(75) 


Since v is any arbitrary function which satisfies Helmholtz's 
equation, v can be replaced by the pressure of an infinite plane 
wave. If this pressure is 

p = Ae 1(k? - wt) 


(76) 
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with e = 0 on a. Equation (77) is the pressure at any point down- 
stream of any aperture upon which plane waves are incident, and will be 
used to describe the diffraction effects of transverse cracks and 
del ami nations upon the sound field produced by the transducer. 

Of course, another (or the same) transducer must be used to sense 
the ultrasonic wave downstream, and the voltage output of this 
transducer is proportional to the integral of the pressure exerted 
over its area, a R . Using the real part convention, the pressure at 
any point at a distance z-| from the sending transducer is given by 


p n = Cz, coswt (coskr + sjnkr )da _ Cz sinwt ( coskr _ sinkr )(Ja 

1 1 J kr 5 r 1 r kr 

a a 

(78) 

The constant of proportionality C is irrelevant since it will be 
eliminated when a ratio is taken later. Equation (78) is rewritten as, 

Pq = A-j coswt - A^sinwt (79) 

in a manner similar to that of Papadakis [10], where A, is times 
the first integral, and A 2 is zv times the second. Thus, the 
proportional response of the receiving transducer is 
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p^ = coswt A^da - sinwt A 2 da 


a R 


(80) 


a R 


Equation (80) is rewritten as 


p.| = coswt - H 2 sinwt 


(81) 


where H-j and H 2 are the first and second integrals, respectively. 
Recalling that the attenuation is obtained from the maximum of p-j , 

Eqn. (81) is differentiated with respect to time, and set equal to 
zero, obtaining 

H-j sinwt + H 2 coswt =0 (82) 


or 


H 2 

tan wt m = ' 


(83) 


where t is the time at which p, is an extremum. From this 
m I 


Si " ^ = (H^H 2 ) 1 ' 2 


(84) 


and 


± H, 


cos wt = 
m 


( H ^ + H ^) 1 '' 2 


(85) 


Substituting Eqns. (84) and (85) back into Eqn. (81), the maximum 
response of the receiving transducer is found to be 


P 1 max = ± (H? + H 2 ) 1/2 


( 86 ) 


The ambiguous sign results from the oscillatory nature of sinusoidals. 
If the procedure is repeated for some other value of z, say z 2> 
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then the attenuation due to diffraction is 

(Pi max) 

“D " 20 109 Tp 2 max) (87) 

where the plus sign has been used in Eqn. (86). 

The actual integrations cannot be evaluated in closed form, hence 
a numerical technique was used. This is discussed next, before 
proceeding to the theoretical and experimental results. 



6. NUMERICAL INTEGRATIONS 

The difficulty with finding solutions to problems in terms of inte- 
grals, such as was done for the diffraction problem, is that these inte- 
grals must eventually be evaluated. This can be done in closed 
form only for very specific cases, and, in general, some other technique 
must be used which, will yield a good approximation of the desired 
integral. It is a general rule that the better the approximation, the 
more complicated is the approximating formula, and vice versa, but a 
reasonable compromise between accuracy and simplicity exists in a 
technique known as Gaussian quadratures. A very good reference on the 
use of quadrature formulae is Kopal [11], and an extension to surface 
and volume integrals is briefly sketched in Zienkiewicz [12]. 

Gaussian quadratures can only be used to evaluate integrals with 
limits of -1 and 1. The approximation for such an integral is 

N 

f (x)dx = l f(a . )H. (88) 

J] 1=1 

where f(x) is the function to be integrated; f(a^) is that function 
evaluated at x - a-; and is the appropriate weighting function 
for a particular value of a.. N is the degree of the quadrature 
formula, and can be any integer larger than two. The H/s and a.'s 
are suitably chosen so that Eqn. (88) is exact when f(x) is any 
polynomial of degree 2N-1 or less. For example, if N were four, 

Eqn. (88) would be exact for polynomials of degree seven or less. 

Values for a^ and H. are tabulated in references [11] and [12] 

N = 2 to 16. ■ 
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Most integrals do not have limits of -1 and 1, but a suitable 
change of variables can transform the integral to one that does. As 
an example, consider the integral 


I = f(x)dx 

x. 


(89) 


V 1 


The linear transformation 


x 0 - x, x 0 + x, 
X = (-2-s — -)n + 


(90) 


will suffice, because as x varies from x-j to x^> n varies from 


-1 to 1. Then dx is given by 


x 9 ~ X-, 

dx = (-^-5 — -)dn 


(91) 


and Eqn. (89) is transformed to 

1 


I = | F ( n) ( 

*■1 


x 2 “ X F 


)dn 


(92) 


where 


x 9 - x, x 0 + x, 
F(n) = f [ {~~~n — ~)n + — -x — -] 


(93) 


In accordance with Eqn. (88) » I is then approximated by 

x 9 - x, N 

I - - ■ -* 1 l H. F(a . ) (94) 

' i=l 7 7 

The change from one-dimensional integration to two-dimensional 
(surface) integration is only slightly more complicated, as long as 
the surface of integration is completely contained within a plane. 
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Consider, first, integrations in polar coordinates. The function to be 
integrated is f(r, e) and the limits on r and e are r-j and 


r 9 , and e, and e«, respectively. Thus 

r 2 0 2 
f f 


• I = 


f(r, e)rdrd0 


r l 0 1 


A change of variables is again necessary, and is performed by 

vy v r 1 + r 9 

r = (-^ T J -)n + -V 1 

0 = (J ’ ) Y + J-2 


(95) 


(96) 


I is then given by 
1 


I = 


r r 

-1 -1 


r ? -.r. r, + r ? r p - r, e ? - 0, 

F(n, y)[( 5 )n + —5 — ) ( o )dndy 


(97) 


Now, let I -j ( y) , a function of y, be given by 


Ii (y) = 


F(n» y)[( 


-1 


V' r i r i + r ? r ? - r i 
— 2 — ~) n + 1 (~^~2 — “) dn 


Then I is given by 


f 0p ~ 0-j 

i - j I] (y)(-=“ 2 — )dy 


-1 

and an application of Eqn. (71) yields 

N 


0 ? - 0. y 

I " (•-%-• L ) l H I,(a.) 
d i=l 1 1 


(98) 


(99) 


( 100 ) 


where is used to represent the order of the approximation (number 
of terms) in y, or e. Butlj(a.) is given by Eqn. (98) when y = a.. 
Hence 
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r 2 - r-j _n 


W l p (a r a i) [ ( *^*2 L ) a j 


r 0 - r, r, + r 0 
2 In- + _U— 2]H, (101) 


2 J j 


whence I becomes 

( 0 2 *" ®] ) (^2 


v, % N N 
r i) „Y n 


I * 


r 2 - r l 


r l + r 2. 


& J, C <V + - J - r - 6 ' 1 

( 102 ) 


The relationship between the terms in Eqn. (97) and Eqn. (102) is 
obvious, upon inspection, and once Eqn. (97) is written, Eqn. (102) 
can be deduced immediately, without need of the intervening steps. 

As a final example, consider the integration of the function 
f(x, y) over the region defined by 


9] U) < y £ 9 2 ( x ) when x ] l x 1 x 2 


(103) 


The desired integral is 


x 2 g 2 (x) 


I = 


f(x, y)dxdy 


x 1 g-|(x) 


(104) 


The variables x and y are transformed by 


x 0 - x, x 9 + X-, 

x = (-^^~ i }n + 


and 


G 9 (n) - G,(n) G ? (n) + G,(n) 

Y = [— 2— 1 lY + — — 2 ~~^ — 


where 
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X 2 X 1 


x 2 + X 1 


6,(11) = g,t( 2 ' ) n + 2 ' ) 


(107) 


From Kreyszig [13] 




(108) 


or 


Xp ~ x, G 9 (n) ~ G,(t)) 
dxdy = (-— — L ) C— 2 — l d ndy 


and Eqn. (104) becomes 
1 1 

I = | | F(n, y) — — — [G 2 (n) - G^rOldndy 

-1 -1 


(x 2 - x^ 


whereupon I is approximated by 
N N 


I 


x 2 - x 


(109) 


(no) 


j- 1 ,!, jJ F < a j’ vw - G i (a j )1H i H j (11,) 


Equation (111), then, is an approximation appropriate for any 
surface integral in the x - y plane. The actual formulas used in 
finding the diffraction losses are given in the appendix, as well as 
the computer programs used to evaluate them. The results of the 
numerical computations are discussed next. 


7. NUMERICAL RESULTS 


To simplify the calculations required to find the diffraction 
losses due to damage development in composites, certain assumptions 
must be made. These assumptions are not justified by any physical 
arguments, but, in as much as this text is meant to be of an exploratory 
rather than of a definitive nature, they are acceptable as a first 
approach to solution. The specimen is assumed to be infinitesimally 
thin, and serves only to cover the transducer in damaged areas 
(locations of cracks or delaminations). Replacing the fused silica 
delay block with an inviscid fluid has already been discussed, but, 
inherent in this substitution is the assumption of no attenuation in 
the quartz, which is only approximately correct. Also implied in 
Eqn . (77) is that the wave is monochromatic, which is, of course, 
not true in pulse-echo work. The next approximation applies only to 
the diffraction effects calculated for transverse cracks, and is 
obtained as follows . 

For high frequency waves, k is much larger than — , except 

r 

for small r, and Eqn. (77) can be simplified to 



izkA 

2tt 


r e i(kr - wt) 


o 



( 112 ) 


Another simplification is afforded if z is assumed large enough so 
that the ratio — is approximately one for the values of r which 
occur in the integrand. Then Eqn. (112) becomes 
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P Q = 


ikA 

2tt 


- wt) 


dS 


(113) 


which is essentially the formula used by Papadakis [10]. 

Since this report is concerned mainly with pulse-echo measurements, 
the sending transducer is assumed to be the same as the receiving 
transducer, though it is displaced in the z -direction by increments 
of 2nL to account for the distance traveled by each echo. But 
more importantly, the receiving transducer is covered by exactly the 
same cracks, in exactly the same location, as the sending transducer. 

Let the area of the transducer be denoted by Ay, and that of 
the cracks by A To differentiate between the sending and receiving 
transducer, a superscript of S or R will be used with the appropri- 
ate area. Then the response of the receiving transducer will be 
proportional to (real part convention) 




,R ft R 


S ,, 5, 


(A^A£) (Ay-A£) 


. i ( kr - wt) 


da S da R 


( 114 ) 


which can be expanded into the form 




e i(kr - wt) 
r 


. S 

da da 



It is a fairly simple matter to show that the second and third 
integrals are exactly equal. The last integral is assumed to be small 
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} 




enough to neglect, and Eqn. (115) is rewritten 


P = 




i (kr - wt) 


da R da S 


e i(kr - wt) R s 
— — — da da 


K Ay 


A S 
T C 


(116) 


Equation (116) was used, instead of the more exact formulation of 
Eqn. (115), to greatly reduce the amount of computing time required. 
The first integral of Eqn. (116) is of special interest because it is 
used to find the diffraction losses from the transducer alone. 


7.1 Accuracy of the Integration Techniques 

The pressure at any point along the axis of the sending transducer 
can be obtained in closed form by integrating Eqn. (77). If z is 
the distance from the transducer to the point on the axis, and if D 
is given by 

D = (R 2 + z 2 ) 1/2 (117) 

where R is the radius of the transducer, then the pressure at any 
point Q (0, 0, z) is 

2 

Pq * [1 + % - 2 | cosk (D - z)] 1/2 (118) 

Because the integration of Eqn. (77) to obtain Eqn. (118) is the 
same type of integration used later to calculate the attenuation due 
to diffraction, this was considered to be an ideal means of checking 
the accuracy of the Gaussian quadrative formulae. This check was 
performed in both rectangular and polar coordinates, since they were 
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to be used in the later calculations. Three values of k (16tt, 20tt 
and 24tt ) were used because they correspond to the range of frequencies 
used in the laboratory (4.8, 6, and 7.2 MHz, respectively). Table IV 
presents the results of the numerical integrations in both coordinate 
systems, as well as the exact answers. The results for all three 
values of k were very similar, and only the results for k = 20tt 
are presented. Note that the polar coordinate formula is very 
accurate for all but the smallest values of z. The rectangular 
formula is less so, but is sufficiently accurate for the values of z 
required (z greater than 1.0 cm). 

7.2 Diffraction by Transverse Cracks 

Because it is necessary to evaluate the diffraction integral 
numerically, rather than in closed form, it is possible to find the 
diffraction losses for only a finite number of cases out of the infinite 
combinations of crack location, length, shape, width, and operating 
frequency which can exist. It was decided, then, to choose a maximum 
of four crack spacings, three crack widths, and three frequencies for 
a total of thirty-six possible combinations, while the cracks were 
assumed to be rectangular in shape, and to extend completely across 
the width of the specimen. The cracks were also assumed to be evenly 
spaced across the width of the transducer so that symmetry could be 
used to reduce the amount of camputation time required. The crack 
spacings used were 0.075 cm, 0.098 cm, 0.141 cm, and 0.254 cm. The 
radius of the transducer used was .635 cm, and for this size transducer, 
the above spacings correspond to 16, 12, 8, and 4 cracks, respectively. 


62 


TABLE IV 

COMPARISON OF THE EXACT FORMULATION AND GAUSSIAN APPROXIMATION 
OF THE PRESSURE ALONG THE AXIS OF THE TRANSDUCER (k = 20-rr) 


z Exact Gaussian Approximations 

(cm. ) Pressure Rectangular Polar 


0 

1.000 

0.000 

0.001 

1.001 

0.116 

0.01 

1.000 

1 .152 

0.02 

0.982 

2.248 

0.05 

0.948 

4.781 

0.10 

1.142 

5.553 

0.20 

1.192 

0.718 

0.50 

0.554 

0.540 

0.80 

1.123 

1.124 

1.0 

0.870 

0.869 

1.5 

1.514 

1.514 

2.0 

0.109 

0.109 

5.0 

1.898 

1.898 

7.5 

1.491 

T . 491 

10.0 

1.182 

1.182 

25.0 

0.501 r 

0.501 

50.0 

0.253 

0.253 

75.0 

0.169 

0.169 

100.0 

0.127 

0.127 


0.000 
0.770 
0.985 
0.982 
0.948 
1.142 
1.192 
0.554 
1 .123 
0.870 
1.514 
0.109 
1 .898 
1.491 
1.182 
0.501 
0.253 
0.169 
0.127 
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across the total width of the transducer. As was mentioned earlier, 
the frequencies used were 4.8, 6, and 7.2 MHz, and the cracks were 
assumed to be 0.008, 0.016, or 0.024 cm. wide. Of the possible 
thirty-six cases, twenty-three were actually used; specifically, 
they were: 

1) 4.8 MHz (k = 16u) 

a) 8 cracks, all three widths 

b) 16 cracks, all three widths 

2) 6.0 MHz (k = 20ir) 

a) 4 cracks, 0.016 and 0.024 cm wide 

b) 8 cracks, all three widths 

c) 12 cracks, all three widths 

d) 16 cracks, all three widths 

3) 7.2 MHz (k = 24*) 

a) 8 cracks, all three widths 

b) 16 cracks, all three widths 

The diffraction curves for the above cases are presented in 
Figs. 14 through 21. The horizontal scale is the distance z traveled 
by the wave, divided by the length of the delay block, 2.55 cm. The 
first echo, then, is at z/L = 2, and the second at z/L = 4. 

Included in each figure is the diffraction loss curve for the trans- 
ducer alone. The following procedure is then used to find the 
attenuation change that would occur for the given damage state. 

The original diffraction loss is found by taking the difference 
in the heights of the transducer loss curve between z/L = 2 and 
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Figure 14. Attenuation vs. z/L for 4.8 MHz, 8 
Cracks 
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Z/L 


Figure 15 . Attenuati on vs. z/L for 4 . 8 MHz, 1 6 
. 'Cracks'. 




Figure 16. Attenuation vs. z/L for 6.0 MHz, 4 Cracks 




ATTENUATION (dB) 
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ATTENUATION (dB) 
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ATTENUATION (dB) 



z/l 

Figure 19. Attenuation vs. z/L for 6.0 MHz, 16 Cracks 
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Z/L 


Figure 20. Attenuation vs. z/L for 7.2 MHz, 8 Cracks 



ATTENUATION (dB) 
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z/L = 4. For example, in Fig. 14, where k = 16ir, this difference 
is 0.21 dB. After 8 cracks, 0.016 cm wide, have completely formed, 
the loss, again found by taking the difference on the appropriate 
curve, is 0.86 dB. Thus, there is a net change of 0.65 dB. The 
results for all twenty-three cases are presented in Tab!*? V. 

Several facts are of note. Firstly, by examining the column 
for the transducer losses, one can see that an increase in operating 
frequency is not always accompanied by a decrease in the initial 
attenuation, as is commonly believed. Secondly, an increase in the 
number of cracks present does not necessarily mean an increase in the 
diffraction losses, as is evidenced by the listings under 6.0 and 
7.2 MHz, although the above is true at 4.8 MHz. And thirdly, 
the net change can be negative; i.e., the apparent attenuation of 
the specimen may decrease as more damage is introduced. 

Even though insufficient cases have been examined to properly 
generalize the above results, Fig. 22 is presented as a graph of the 
expected attenuation change versus crack spacing with frequency and 
crack width as parameters. No strict conclusions should be drawn 
from this graph until more data has been obtained. It is shown here 
merely as an aid to visualization of the results in Table V. 

From these results, it is readily apparent f^t frequency plays 
a major role in the amount of change which can be expected. By 
simply varying the frequency from 4.8 to 7.2 MHz, the diffraction 
loss for one case would change from 1.74 dB to -0.14 dB, a net 
difference of almost 2 dB (Recall, though, that the bond loss term 
can also vary quite a bit with frequency, overwhelming more minor 
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TABLE V 

DIFFRACTION LOSSES OF VARIOUS CRACK SPACINGS AND WIDTHS 


Freq. 

(MHz) 


4.8 


6.0 


No. of 
Cracks 

Crack 

Width 

(cm) 

Loss from 
Transducer 
(<JB) 

' Loss from 
Cracks 
(dB) 

Net 
Loss 
... (dB) 

8 

.008 

0.21 

0.86 

0.65 


.016 


1 .39 

1.18 


.024 


1.93 

1.72 

16 

.008 


0.93 

0.72 


.016 


1.51 

1.30 


.024 


1.95 

1.74 

4 

.016 

0.20 

0.91 

0.71 


.024 


1.25 

1.05 

8 

.008 


0.54 

0.34 


.016 


0.81 

0.61 


.024 


1.04 

0.84 

12 

.008 


0.46 

0.26 


.016 


0.63 

0.43 


.024 


0.73 

0.53 

16 

.008 


0.40 

0.20 


.016 


0.49 

0.29 


.024 


0.49 

0.29 

8 

.008 

0.59 

0.71 

0.12 


.016 


0.79 

0.20 


.024 


0.83 

0.24 

16 

.008 


0.59 

0.0 


.016 


0.54 

-0.05 


024 


0.45 

-0.14 


75 


changes in the diffraction loss). Obviously, the frequency at which 
these tests are conducted in the future should be strictly controlled 
and accurately measured. 

7.3 Diffraction by Delaminations 

To calculate the diffraction effects of delaminations, it was 
assumed that a transducer with a diameter of 0.5 in. was centered 
on a specimen 1 in wide. The boundaries of the delaminations 
were assumed to be straight lines, parallel to the edges of the 
specimen, which have progressed to equal depths on either side of 
the specimen. Of course, until the delaminations are a quarter of an 
inch deep, they will have no effect on the transducer, and the apparent 
attenuation due to diffraction is the same as for the transducer 
itself. If the delaminations are a half of an inch deep, they 
completely cover the specimen, and no sound will pass through. 

The diffraction loss curves for various depths of delaminations 
are shown in Fig. 23 for k = 20ir. The expected losses for the 
1 inch long delay block were calculated as detailed above, and are 
presented in Fig. 24. Note that the expected losses are fairly 
insignificant until the delaminations penetrate to about 0.29 in. , 
but that they rise fairly quickly after that point, reaching 2.8 dB 
when the depth is 0.4 in. The predicted loss is around one-half of 
a decibel when the depth of penetration is only 0.31 in., again 
attesting to the speed with which the attenuation rises once the 
delaminations penetrate beyond the critical depth of 0.29 In. 

The fact that delaminations have not generally been found to 









rr 





CHANGE 



I DEPTH OF DELAMINATION (IN) 

Figure 24. Change in Attenuation vs. Del ami nation Depth 




penetrate more than a quarter of an inch into the width is not a 
serious problem, since the transducer is usually not well centered 
on the specimen. A displacement of a tenth of an inch from the 
center is not unusual (in fact, it is quite likely) and would result 
in extensive coverage of one side of the transducer. 


8. CONCLUSIONS AND SUMMARY 


From this study, it is apparent that the existence of transverse 
cracks or delaminations can indeed account for the changes in 
attenuation which occur as specimens are loaded in quasi-static 
tension, and thus, these changes could serve as an indicator of 
damage formation and growth in the specimens in regions not readily 
accessible to visual observation. Moreover, the predicted changes 
are well within the range of values found experimentally, and also 
show both increasing and decreasing trends evidenced by experiments. 
Both the amount end the direction of the change are closely related 
to the frequency used in the test. In the future, the frequency 
should be closely controlled. 

It is the as yet unproven contention of the author that the 
very sudden rises in attenuation are caused by damage formation, 
such as a transverse crack appearing, while the more slow and 
gradual changes are due to such phenomena as the already existing 
cracks opening wider under load; and that both the rapidity and the 
amount of such changes can give some details as to damage formation 
and growth. 

Because of the general nature of this work, little emphasis was 
placed on obtaining exact correlations for specific specimens. 

Studies of this type are certainly in order, as are also investigations 
of errors, such as improperly functioning transducers, which can 
mask or distort the experimental data. 

With the improvements which are sure to arise from additional 
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investigations, the method described here-in appears to be of some 
aid to investigators studying the formation and growth of damage in 
thin composite laminates. 
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A. APPENDIX 


The integral equations described earlier, along with their 
approximation formulae are listed below. 

A.l The Pressure Along the Axis of the Transducer 

The maximum pressure at some point lying on the axis of the 
sending transducer is 


p = (H 2 + h |) 1/2 


In radial coordinates, and H 2 are given by 

R 


ti _ f /coskr , sinkn j 

1 ( 17~ "T' p p 

0 xr r 


and 


H 2 = zk 


, coskr sinkrv . 

( 5 ” o )pdp 


kr‘ 


The approximations for the above are 


zkR 1- u rcoskr , s inky's 

i ~ -T~ l h. (1 + a.) ,[ 3- + —o“J 

1 4 i=1 1 1 kr d r 


2 16 


coskr . sink 


and 


zkR' 


2 16 


- T- l h. (1 + a.) 

4 i v v L 2 


coskr sinkr 


■] 


i=l 


kr' 


where a. and h.. are the Gaussian quadrature constants. The 
numerical value of r is 
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(119) 


( 120 ) 


( 121 ) 


( 122 ) 


(123) 
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r - [f (1 ♦ a*) ♦ W 2 


(124) 


In rectangular coordinates, and are approximated by 


H, - 


7kR 2 16 16 

l l h. h. (1 

2,r i=i j=i 1 j 


„2J/2 r coskr , sinkr-, 

V [~~3“ 2“ j 

kr r 


and 


0 25) 


Ho - 


zkR 


2 16 16 


9 l l h. h. (1 - at) 1 ^ 
2. !=■) & i 3 r 


2,1/2 r coskr sinkr 


kr' 


•] 


(126) 


where 


r R 2 2 . R 2 ,, .2, 2 . 2 -,1/2 

r = [ r a i + r (1 - a i )a j + z ] 


(127) 


A. 2 The Pressure on the Receiving Transducer with No Damage 

Using the approximation detailed in Section 7, the average 
pressure on the receiving transducer is 


P * (H 2 + H 2 ) 172 


(128) 


where 


H, = 2tt 


2-ir R R 

f coskr 


pdpridr|d0 


(129) 


and 


0 0 0 


2ir R R 


H 0 - 2tt 



sirikr 


pdpndnds 


(130) 


0 0 0 


where 


r = [(ncose - p)^ + n 2 sin 2 e + z 2 ]^ 2 (131) 


The approximation formulae are 


d 4 2 16 2 16 16 

R ir 


coskr 


H, - Sjf- I I l l 0 + a-)(l + a m ) h. h h aa 
' 16 1»1 j=l riFl n=l 1 m .1 m n r 


(132) 


and 


0 4 2 16 2 16 16 

R 7T 


sinkr 


9 “ I I I Z 0 + a . ) (1 + a ) h. h h ^ 

2 16 i=l j=l n=l 1 m 1 m n r 


(133) 


where 


r = {£- [(1 + a m ) cose jn 


(i +a i >] 2 + r (i +a / A 

(134) 


and 


e jn = J[(2j-l) + a„] 


(135) 


A. 3 Response of the Receiving Transducer in the 
Presence of Transverse Cracks 

As discussed in Section 7, H-j and were modified by 

9 ■■ , l. 

subtracting the integrals denoted by and H ^ below. Otherwise 
Eqn. (128) still applies in finding the average pressure on the 

I . 1 

receiving transducer. H-j and are 




Kjoasea 
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m\ 
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2 tt R 


H, = 2 



pd pd8da C 


(135) 


A c ° 0 


2 tt R 

H 2 = 2 M f ^ > d » ded ° C 

L o o 

V/. 


(136) 


The integration over the area of the cracks was done with a 16-point 
quadrature formula in the y direction and a 4-point formula in the 
x direction. Let x. denote the centerline of a crack, and w 
the width of a crack. The height of the crack is y. , and is 
given by 


y f = (R 2 - x 2 ) 1 / 2 


(137) 


The approximations for H-j and LL, are 


H 


d 2 NC 16 4 16 16 2 

ttk w 


coskr 


i « TT* l l 1 l l y h. h h h (1 + a )y. ^ 

1 8 1=1 j=l <=] m-_l n-1 s=l J 1 ™ 1 r 

(138) 
sinkr 


^2 NC 16 4 16 16 2 


V T* T l l l l I H. h t h h o + a )y 
° i=l j=l £=1 ffi=l n=l s=l J £ m n m 1 


where NC is the number of cracks. To formulate r, let 


(139) 


R 1 = 2 (1 + a m ) cose s[J 


(140) 


and 


r - & 
K 2 " 2 


+ 3 m> sin0 sn 


041) 


.1 ' 




C\».S 'af' :j. v rz a# T ■ 


1 2 it 


r - [Rj + (R 2 - a, y 1 ) 2 + z 2 ] 1/2 


(143) 


e sn = | [2. - 1) + a n ] 


(144) 


A. 4 Response of the Receiving Transducer in 
the Presence of Delaminations 


The response of the receiving transducer is also governed by 
Eqn. (128). The formula for and H 2 in rectangular coordinates 


□ _ zk 

H 1 2n 


R 1 (R 2 -n 2 ) 1/2 R 1 (R 2 -x 2 ) 1/2 

[I [ (5M£. + £2skC) dndydxdy 

\ J (R 2 - n 2 ) 1/2 -R, i(R 2 -x 2 ) 1/2 ^ kl " 

(145) 


u - zk 

H 2 2tt 


R 1 (R 2 -n 2 ) ,/2 R 1 (R 2 -x 2 ) ,/2 


rcoskr sinkn 


in, V-nV /2 \ -(R 2 -x 2 )^ 2 


( 146 ) 


where R^ is the perpendicular distance from the center of the 
transducer to the inside edge of the delamination, and 


r/ >2' , x2 A 2 n l/2 

= C(x - n) + (y - y) + z ] 


88 


and 


The approximations for H-j and are 


zR,k 16 16 16 16 


H i = -r- i £ £ £ ft, h . h h D D. (^£ + 

1 2 * i=l j=l H=1 iifI 1 : * "I 1 J r 2 kr 3 


( 148 ) 


zR 2 k 16 16 16 16 k 

H, = — T T T 1 h. h. h h D. D. (£2|£E 
2 2ir 1 i, ji, 1 J 4 n, , j 2 . 


where 


sinkr> 

— 

kr 


(149) 


D i = (R 2 - R 2 a 2 ) 1/2 (160) 

and 

r = [R, (aj - a ( ) 2 + (D^ a s - D. a k ) 2 + z 2 ] 1/2 (151) 

A. 5 List of Symbols used in the Programs 

The following list defines the relationship between the symbols 
in the computer programs and the nomenclature in Eqns. (119) through 
(151) 

A (I) : a. (16-point Gaussian quadrature) 

H (I) : h. (16-point Gaussian quadrature) 

Z (I) : z 

BETA : k 

R : R 

HI : H 1 
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H2 

: H 2 

X (I) - 

: x i 

Y (I) 

: *i 

ABAR (I) 

: ' a 1 

HBAR (I) 

: h, 

HIT 

: h; 

H22 

: H 2 

P 

: P 

RHO 

: r 


(4-point Gaussian quadrature) 
(4-point Gaussian quadrature) 


A. 6 Computer Programs 


The following is a listing of the computer programs used in 
computing the appropriate values of attenuation. 


o o 


C 

C 

C THIS PROGRAM CALCULATES THE PRESSURE ALONG THE AXIS Cf THE 


C TRANSDUCER IN X-Y COORDINATES, 

C INPUT DATA: 

C H(U) - HI 16) (GAUSSIAN QUADRATURE) Fib. IS 

C A IU J - A l 16) (GAUSSIAN QUADRATURE) FIS. ip 

C BETA (WAVE NUMBER/ PI > FIB. 15 

C R (RADIUS OF THE TR ANSDUCtR ) FIS. 15 

C NZ (SIZE OF THE ARRAY Z) 12 

C Z( l) - Z ( N Z ) FIS. 15 

C 

C NOTE: ALL DI HEW SIGNS ARE IN CM. 


IMPLICIT REALMS (A-H,0-Z> 

DIMENSION At 16) ,H(16> , Z( 50) 
PI=4.D0*DATAN(1.D0) 

READ ( 5*500.) H»A,BETA*R 

READ( 5,510) NZ 

READ (5, 500) { Z ( I ) , 1=1, NZ) 

BETA=BETA*PI 

500 FORMAT ( F 18 . 1 5 ) 

510 FORMAT ( 1 2 ) 

WRITE (6,600) BETA 

600 FORMAT ( • i* , 5X, • BETA = • , F12 . 6//// / / ) 
DO 20 1=1, NZ 
H 1=0 • DO 
H2=Q . DO 
DO 10 J=1 ,16 

B=DSQRT(1,D0— A(J) *#2 ) *H( J) 

R1=R*R*A( J)**2 



i)Q 10 K=1 » 1 6 

RHu=DSQRT (kl-“ ( 1 . DO-Al K) *~2) +R*R*A ( K)**2+Z 1 I l 1 —/. ) 

R2=RH0*RHG 
R3=RHO*R2 
C=OCOS( BETA^RHO } 

S=OSIN( BETA*R HO ) 

H 1= H 1 +8 * H ( K ) * ( C /BETA / R3 + S / R 2 ) 

10 H2=H2+3*H(K)*(C/R2~S/rtcT4/R3) 

H1~H1*ZCI )*R*R/2.DD/PI*BtTA 
H2=H2*Z(I ) ^R*R/ 2»D0/PI*BcTA 
P=DSQRT(Hl-*2+H2**2 > 

Q=DSORT(R*R+Z( I )**2) 

Pl*DSQRT( l.DO+ZI I )**2/Q/it-Z.00*Zl I ) /D*DCUS ( BE TA* ( 0- Z( I) ) ) ) 
610 FORMAT (5X» 3< F10 «6»7X) > 

20 WRITE ( 6 »6l j ) Z ( I ) , HI , H2 t P , P 1 
STOP 
END 



O U U U O o o 


THIS PROGRAM CALCULATES THE PRESSURE ALONG THE AXIS PT THl 
TRANSDUCER IN RADIAL COORDINATES. INPUT DATA IS THE SAME • 
AS IN THE PREVIOUS PROGRAM. 


IMPLICIT REALMS {A-H,0-i> 

DIMENSION A (16) ,H(16) ,Z(50) 
PI=4.DQ*DATAN(1.DQ) 

READ l 5 1 500 ) H f A t BET A, R 
READ ( 5 , 51 0 ) N l 
BET A=8ETA*P I 

READ 15,500)IZII),I=I,NZ) 

500 FORMAT (FI 3* 15) 

510 FORMAT ( 12 ) 

00 5 1=1,16 
5 A( 1 1 — A ( Ii + l.DO 
WRITE ( 6, 600 ) 

6 00 FORMAT ( • I » > 

WRITEt 6, bZQ) BETA 

620 F0RMAT(5X, • BETA = • ,F12 . 6///// ) 

DO 20 1=1, NZ 
H1=0 .DO 
H2=0 . DO 
DO 10 J=1 ,16 

RHO==DSQRT (Z( I)**2+A( J )**?*R*R/4. 00 > 

R2=RHQ*RH0 

R3-R 2*RH0 

C=DCOS(BETA*RHOi 

S=DSIN<8£TA*RH0) 

H1~HL+A(J )*( S/R2+C/BETA/R3)*H( J) 


10 H2-H2+A ( J )*( C/R2-S/B=TA/R3> *Hl J) 

Hl=Hl-SfcTA*R/2.Q0*Z( I )*R/2.D0 
H2=H2tbETA-R/2.O0^2i I)*R/2.D0 
P=DSQRT(H1**2+H2**2) 

O^DSQRT ( R *R+Z ( I ) =**'2 ) 

P1=DSQRT ( 1 ,00+Z ( I )*^2/D/0-2.D0^Z( I ) *DCO$ l BETA* t 0-2 ( I) ) ) /O) 
610 FORMAT l5X f 5(F10.6,7X) ) 

20 WRIT216,6L0) Z( I) »H1,H2,P,P1 
STOP 
END 



c 

c 

C THIS PRQGRAi'1 CALCULATES THE APPARENT ATTENUATION AT 
C VARIOUS DISTANCES, WITH CRACKS PRESENT. 

C SUBROUTINE SU31 MUST BE USED IN CONJUNCTION WITH THIS 


C PROGRAM. 

C INPUT DATA: 

C Hll) - H ( 16 ) (GAUSSIAN. QUADRATURE) FLB.L5 

C A ( 1 ) - A ( 16 ) (GAUSSIAN QUADRATURE) FIB. 15 

C R (RADIUS OF THE TRANSDUCER) FIB. 15 

C W (ANGULAR FREQUENCY) FL8.15 

C V (VELOCITY OF SOUND! F18.15 

C WID (WIDTH OF THE CRACKS) F18.15 

C NC (HALF THE NUMBER OF CRACKS) 12 

C H3AR ( 1 ) - HEAR (4) (GAUSSIAN QUAD.) FIS. 15 

C A8AR ( 1 ) - ABAR14) (GAUSSIAN QUAD.) Fib. 15 

C NZ (SIZE UF THE ARRAY Z) 12 

C Z( 1) - Z(NZ) FIB. 15 

C 

C NOTE: ALL DIMENSIONS ARE IN CM. 

C 

C 

IMPLICIT REALMS (A-H,0-Z) 


COMMON E, HI, H2, BETA, H( 16) , A( 16) ,P, S(2, 16) ,C(2, 16) , PI 
COMMON S2 (2,16) ,C2(2»16) 

DIMENSION X(50) ,Y(50) , Z ( 100 ) , P 1 ( 100) , P2 ( 1 00) , HEAR ( 4 ) , ABA* ( 4 ) 
READ (5,500) H , A ,R , M , V , W ID 
530 FORMAT ( F18. 15 ) 

RE AD( 5 » 510 ) NC 
510 FORMAT (12) 

RE A0( 5 , 500) HBAR,A3AR 
Pl=4. DO^DATAN ( L.DO) 


B ETA- w’/ V 
DtLX=R/(NC+.5> 

X< L)=DELX/2.D0 
X ( 1+NC I ~ — X ( 1 ) 

Yil)=D$QRT(R*R-X( 1)**2> 

Y ( l + NCI —Y t L J 

DO 5 I-2»NC 

X'{ I) = X( I-1)*QELX 

Y( I )-QSQRTlR*R-X(I )**2) 

X ( I +NC ) =— X( I) 

5 Y ( I+MC)-Y ( I ) 

NC2=2*NC 

WRI Tfc ( 6 » 61 ) 

61 FORMAK • L» ,TL4, »Z* ,T24,' S' ,T43,«H1* ,T73,'H2'/) 
READt 5,510) NZ 
READ (5, 500) (Zd ), 1 = 1, NZ) 

DO 6 K=l, 2 
T2=P I /4. DO 
T i=P I /A.DO* ( 2*K— 1 ) 

T3=T1/2.D0 
T 4=T 2/2 .DO 
DO 6 N-1,16 
TNl=Ti+A(N)*T2 
TN=T3+A(N)*T4 
S2( K, N)=DSI Nt TNI) 

C21K, N)=DC0S(TN1) 

S ( K , N ) = D S I N { T N ) 

6 C(K,N)-DCOS(TN) 

DO 30 L=l, NZ 
CALL SUBHZ(L)) 

P11L)=P 

HI 1=0. DO 



cn 


>122= a. DO 
DO 20 K= 1 » 2 
DO 20 N=1 , 16 
DO 20 H= 1 » 1 6 

Rl=R/2.00*( l.DiH-A(M) ) *C(K,N) 

R2=R/Z.00*< 1.00+AfM) )*S(K,N> 
u=H( M)*H(N) *( 1 .DO+A(M) ) 

DO 25 11=1,4 

k 3=R L-X ( I >-WID/2.D0*A6ARt ID 
61 = S*Y( D*HBAR( III 
DO 25 1 = 1 , i'4C2 
DO 25 J=1 ,16 

RHQ=DS QRT (R3**2+Z { L ) **2 + < R2-Y C II^A ( J ) i **"2 ) 

C i=DC OS ( 6ETA*RHG) /RHO 
S 1=DS IN( BET A^RHQ ) /RHO 
HH=H11+HC J)*B1*C1 
25 H22=H22+H( J ) *B1*S1 
20 CONTINUE 

H 1 =H1 — P I * R* R= W I 0/ 8 . DO* HI 1 
H2=H2-PT*R«R^W1D/8.D0*H22 
WRITE (6»t>2) HI, H2 

62 FORMAT (T30,D23. 16, T63,D23 .16/ ) 

30 P2( L)=DSQRT( Hl**2+H2**2) 

WRIT £(6,63) 

63 FORMATS * 1' , T14, * Z* »T25» ' ALPHA 1 ' , T45 , * ALP HA 2 * , T70, * PI * ,T 50 , ' P2 * / ) 
A 1=P I *R*R 

PB=A1/6ETA*2.D0*PI 
PA=AD2.D0*PI/8ETA 
DO 40 1=1, NZ 

AL1=20 • DO*DLOGlO ( PB/PI ( I ) ) 

AL2=20.D0*0LGG1 0( PA/P2 ( I ) ) 

64 FORMAT (T 10, F8 . 5 ,T2 5 , F 12 . 8 , T 45 , F 12 . 8 ,T65 , F 12.8 ,TB5,F12.8) 


CM* WIO 5 -'.*) 


40 WRI TE ( 6 ,64) Z ( I ) , AL 1 , AL2 , PI II ) , P2 ( I ) 

WRITE ( 6,70) NC2 , WI L) T 

70 F0RMAT(///T15, 'THERE ARE • ,13,* CRACKS ,*» FO .6 , • 

STOP 
END 

SUBROUTINE SUBl(Z) 

IMPLICIT RHAL*B{A-H,0-Z) 

COMMON R, Hi , H2 , BETA » HI 16) , A (16) »P,S(2,16) ,C ( 2 , 16) ,PI 
COMMON S2C2,16) ,C2( 2, 16) 

00 5 1=1,16 
5 A(I )=A( n+i.oo 
10 HI =0 • DO 
H2=Q .00 
Z2=Z**2 

S1=Z*2.D0*PI/6ETA/R**2 

DO 20 K=L,2 vo 

DO 20 M=1 ,16 ^ 

DO 20 N= 1 , 1 6 
RH01= AIM) *C2IK,N) 

RH02=( A<M>*$2(K, N)*R/2.D0 )**2+Z2 
8=A(M)*H(M)*HIN) 

DO 20 1=1,16 

RHO=DSQRTIR**2/4.DO*(RH01-A( m**2+RH02) 

H1=H1+A{ I )*H{ n*B*DCaS{ B£TA*RHO ) /RHO 
20 H2=H2+AU)^Hl I ) *3*DS IN I BETA*RHO ) / RHQ * ' 

Hl=Hl*R**4*Pl*PI /16.D0 
H2=H2*R**4^P I *P 1/16. DO 
P=DSQRT(H1**2+H2**2 ) 

WRITE i 6, 61) Z , SI , Hi , H2 

61 FORMATTTiO, F8.5 ,T20»F8. 5,T3Q*D23. 16 ,T60,D 23. 16) 

DO 30 1=1,16 
30 AC I ) = A( I ) — l.DO 
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RETURN 

E-MO 


VO 

CO 


c. 

C THIS PROGRAM CALCULATES TnE ATTENUATION IN THE PRESENCE 


C OF DEL AMI NAT I ONS . 

C INPUT DATA: 

C H ( 1 ) - H( 16) (GAUSSIAN QUAD.) Flo. 15 

C A ( 1 ) - A(16) (GAUSSIAN QUAD.) F13.15 

C BETA (WAVE NUMBER/PI) FI8.15 

C R (RADIUS OF THE TRANSDUCER.) FIS. 15 

C R1 (DISTANCE FROM CENTER OF TRANSDUCER 

C TO THE EDGE OF THE DELAMINATION) FI 8.15 

C .\Z (THE SIZE OF THE ARRAY Z) 12 

C Z ( 1 ) - Z(NZ) FIB. 15 

C 

C NOTE: ALL DIMENSIONS ARE IN CM. 


IMPLICIT REAL *8 (A-H.O-Z) 

DIMENSION AU6I *H( 16) ,Dt 16) ,Z(25J 

PI=4.D0*DATAN(1.D0) ' 

READ ( 5 * 500) H*. A* BETA* R, R l 

READ( 5 , 51 0) NZ 

READ (5*500) ( Z ( I ) » 1= i * NZ ) 

50 0 FORMAT (FI 8. 15) .* v - 

510 FORMAT (12) " r 

BET A=B ET A*P I 
DO 5 1=1,16 

5 D( I ) = OSQRT ( R^R-Rl^R 1*A( I)**2) 

WRITE (6*600) BETA* R l » R • ' ‘ 

600 FORMAT* ,5X,' BETA = • ,F12 .6//5X , • Rl = * » F12 .6, 5X * * R = ',F12.6//) 
DO 20 M=1 » NZ 
HI=0. 00 


I 


H2=0 . DO 
Z1=Z1H) 

DO 10 1=1,16 
b=H{ I)*D( I ) 

DO 10 K=1 ,16 
E=B*H( K)*D{ K.) 

RR=R1*R 1* l A (K )-A ( I ) } **2+Z 1*Z 1 
DO 10 J=1 ,16 
F=S*H ( J ) 

DO 10 L=l,16 

RH0=DSQRT(RR+(0(K)*A(L)-D(I )*A( J> ) **2 ) 

R2=RH0*RH0 

R3=R2*RHQ 

C=OCOS t BETA^RHO ) 

S=DS I N( BETA^RHO ) 

H1=H1 + F^H( L ) HS/R2+C/BETA/R3) 

10 H2=H2 + F*H(L )^(C /R2— S/ BETA/ P.3) 
HI=H1*R1*R1-Z1*8ETA/2.D0/PI 
H2=H2*Rl*Rl*Zi-8ETA/2. DO/PI 
P=DSQRTIH1**2*H2**2 ) 
ALPHA=20.D(3*DLOG10(PI*R*R/P) 

610 F0RMAT(5X,5(F12.6,7X) ) 

20 WR I Tc (6,610) Zl, HI, H2,P, ALPHA 
STOP 
END 


